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Unsteady Hybrid Vortex Technique for Transonic Vortex Flows 
and Flutter Applications 

Osama A. Kandil 

This report covers the progress of research work performed under this 
grant in the period of February 16, 1986 to October 15, 1986. During this 
period of eight months, the following tasks have been accomplished: 

1. Integral Equation Method of Full Potential Equations 

The integral equation method has been extended to treat transonic flows 
around airfoils. This is achieved by adding a, volume integral term 

(corresponding to the full compressibility effect) to the classical 
surface integral solution. The gradient terms in the volume integral term 
are calculated by using mixed type differencing which is consistent with 
the mixed nature of the transonic flow. The solution is obtained through 
successive iteration cycles for subcritical flows. In addition and for 
transonic flows, mixed type differencing is used and the Rankine-Hugoniot 
relations are satisfied across the captured shock. Moreover, shock panels 
are explicitly introduced and fitted to sharpen the captured shock. 

The method is applied to airfoils in subcritical and critical flows and 
the results are in good agreement with the experimental data and the 
finite-difference solutions of the full-potential and Euler equations. 
The method has been applied eariler to transonic vortex-dominated three- 
dimensional flows using a low-order distribution of vortex panels. 
Currently the method is extended to solve for unsteady transonic flows 
around pitching airfoils and to solve for the steady and unsteady 
transonic flows around three-dimensional wings. 

The following papers have been published or submitted for presentations 
( copies are attached to this report); 

a. "Transonic Flow Computation Using the Integral Solution of the Full 
Potential Equations," Kandil, 0. A. and Hong, H., submitted for 
presentation at the AIAA Fluid and Plasma Dynamics Conference, 
Honolulu, Hawaii, June 8-10, 1987. 
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b. "Finite-Volume and Integral -Equation Techniques for Transonic and 
Supersonic Vortex-Dominated Flows," Kandil, 0. A., Chuang, A. and 
Chu, L-C., ICAS-86-1,5.4, London, England, September 7-12, 1986, pp. 
629-640 (Vol. 1). 

c. "Transonic Vortex Flow Past Delta Wings-Integral Equation Approach," 
AIAA Journal, Vol. 24, No. 11, 1986, pp. 1729-1736. 

2. Euler Equations Solution Using Central -Differencing Finite Volume Solver 
with Second and F^ourth Order Dissipation terms 

During this period, we have implemented Jameson's Central -Differencing 
Finite Volume Method, which uses Runge Kutta time stepping with second and 
fourth order explicit dissipation terms, into two computer codes tosdve 
for vortex dominated flows over a wide range of Mach numbers which 
includes subsonic, transonic and supersonic regimes. 

The first computer code solves for the supersonic vortex- dominated flows 
which is an exact conical flow for steady inviscid problems. The unsteady 
Euler equations are used to study conical supersonic vortex -dominated 

flows around delta wings with sharp and round edges. For sharp edges, 

separated flow solutions have always been obtained while for round edges, 
separated and attached flow solutions have been obtained. The solution 
for round edges is dependent upon the level of numerical dissipations and 
the grid fineness. We are the first to conclusively show the solution 

dependence upon the numerical dissipation for round edges. Recently, we 
have shown two additional parameters upon which the solution depends. 
These are the total pressure loss due to the entropy production and the 
accuracy of enforcing the boundary conditions. 

The second computer code solves for the three-dimensional flows in 
subsonic, transonic and supersonic flows. The code is fully vectorized 
and is very efficient (For 130,000 grid points, it takes 3,000 CRU for 900 
iteration cycles until convergence — average residuals of 10"'^ - on the 
VPS 32). The code has been used to solve for transonic and subsonic 
vortex - dominated flows and for supersonic vortex- dominated flows as 
well. The results compare well with the experimental results at the chord 
stations downstream of the 18 % chord station with 130,000 grid points. 
Fine grids (300,000 points) are being tested. 
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Work Is underway to solve for the unsteady pitching and rolling 

oscillations with time accurate schemes. 

The results of this work have been reported in the following papers 

(copies are attached with this report); 

a. "Numerical Dissipation Effects in Finite-Volume Euler Solutions for 
Conical Vortex-Dominated Flows," Kandil, 0. A. and Chuang, A., 
Springer Verlag, Computational Mechanics, Theory and Applications, 
Edited by G. lagawa and S. Atluri, Tokyo, Japan, May 1986, pp. XI57- 
XI. 64 (Vol. II). 

b. "Influence of Numerical Dissipation in Computing Supersonic Vortex- 
Dominated Flows," Kandil, 0. A. and Chuang, A., AIAA 86-1073, May 
1986. Accepted for publication in AIAA Journal. 

c. "State of the Art of Computational Vortex Flows," Kandil, 0. A. and 
Chuang, A., SIAM National Meeting, Boston, Mass., July 1986. 

d. "Finite-Volume Euler and Navier-Stokes Solvers for Three-Dimensional 
and Conical Vortex Flows Over Delta Wings," Kandil, 0. A. an Chuang, 
A., AIAA 87-041, January 1987. 

3. Presentations and Briefings : 

The P.I. has given three presentations and briefings at the UAB-NASA 

Langley Research Center; 

a. Briefing to Dr. Eligh Turner, Wright Patterson AFB, on the Vortex- 
Dominated Flow work. 

b. Briefing during Division and Directorate Annual Review of the UAB 
Research work. 

c. Presentation to McDonald Douglas visitors to NASA Langley Research 
Center on the Steady and Unsteady Vortex- Dominated Flow research work, 
September 17,1986. 
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4. Organizing and Chairing Sessions on Vortex Dominated Flows 

The P.I. has organized and chaired the following sessions in the vortex 
flow area; 

a. "Vortex Flows," AIAA 24th Aerospace Sciences Meeting, January 1986. 

b. "Computational Vortex Flows," SIAM National Meeting, July 1986. 

c. "Vortex Dominated Flows I & II," AIAA 25th Aerospace Sciences Meeting, 
January 1987. 
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Extended Abstract 


Transonic Flow Computation Using the Integral Solution^ 
of the Full Potential Equation 


Osama A. Kandil* and Hong Hu 
Department of Mechanical Engineering and Mechanics 
Old Dominion University, Norfolk, VA 23508 


Abstract 


The well developed surface panel method is extended by adding a volume 
integral term allowing the calculation of the full effect of 
compressibility. The full effect of compressibility is calculated by using 
mixed type finite difference scheme consistent with the mixed nature of 
transonic flow. The solution is obtained through successive iteration cycles 
for subcritical flows and for critical flows the solution is obtained through 
satisfying the Rankine-Hugoniot relations across the captured shock in 
addition to the successive iteration cycles. Shock panels are introduced to 
sharpen the captured shock. The method is applied to airfoils in subcritical 
and critical flows and the results are in good agreement with the experimental 
data and finite-difference solutions of the full -potential and Euler 
equations. 


t This research work is supported by NASA Langley under Grant No. NAG-1-S91. 
* Professor, AIAA member. 

Ph.D Graduate Research Assistant, AIAA member. 


1. INTRODUCTION AHD BACKGROUND 


At the present time, computational inviscid flow prediction methods 
available to the aerodynamic designer can be devided into two types: first, 
the integral equation methods more commonly known as panel methods or field 
panel methods of which the first was that developed by Hess and Smith [1] and 
second, field methods (which encompass Finite-Difference, Finite-Volume and 
Finite-Element methods) which have been extensively developed for transonic 
flow following the work of Murman and Cole [2]. 

A brief review of the state of the art of computational transonic 
inviscid flows will be presented here. 

(a) Field Methods in Steady Inviscid Transonics 

The research for inviscid transonics has been made in major three ways: 
(i) the transonic small -disturbance (TSO) formulation, (ii) full -potential 
(FP) formulation and (iii) Euler equation formulation which includes the 
rotational ity effect. 

The field methods of transonics are the new ones, having essentially 
started in 1970. Murman and Cole [2] were the first to achieve a stable 
transonic solution for the two-dimensional TSO equation by using the concept 
of type-dependent differencing. Soon after, this procedure was extended to 
three dimensional flows for swept-wing calculations by Ballhaus and Bailey [3] 
and wing-cylinder calculations by Bailey and Ballhaus [4]. At about this same 
time, numerical procedures for solving the transonic FP equation were being 
developed with suitable mapping procedures. Notable contributions are due to 
Stegar and Lomax [5] and Garabedian and Korn [6], both with nonconservative FP 
formulations for airfoil configurations. The first three-dimensional 
calculation (nonconservative) was introduced by Jameson [7] and was used to 



solve the transonic flow about wings. Subsequently, a solution procedure for 
the conservative form of the FP equation was introduced by Jameson [8] and 
extended to three dimensions by Jameson and Caughey [9, 10], 
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Much of the recent work in steady transonic potential field methods has 
been devoted to finding better ways of solving the full potential equation in 
conservation form and to the development of more computationally efficient 
relaxation schemes, such as fully implicit approximate factorization (AF) or 
multigrid. Recently, several methods which use the strong conservative form 
of the unsteady Euler equation have been developed to solve steady and 
unsteady transonic flows. A typical method for steady, two-dimensional 
transonic flow is developed by Jameson [11]. Basic research involving the 
Euler formulation for steady transonic applications is currently increasing 
and is expected to continue to do so as the potential algorithms enter a more 
production-oriented research phase. 

(b) Integral Equation Methods in Steady Inviscid Transonics 

Although a great deal of progress has been made in solving nonlinear 
fluid flow problems by finite-difference methods, these methods have not yet 
proved to be easily adaptable to complex three-dimensional surfaces. The 
major technical obstacle to computing inviscid transonic flows about complete 
aircraft is the difficulty in generating suitable grids. More recently, the 
difficulties and relatively larger computation time associated with field 
methods prompted a number of workers to reconsider the application of integral 
equation methods to transonic flows. 

Panel methods for subsonic and supersonic aerodynamic applications have 
been in use in the aerospace community since the 1960s and have become 
indispensable tools in aerodynamic analysis and design. 



In the low-speed regime, methods which are directly or indirectly 
obtained from the Green's function solution (lEM) have been developed for 
steady and unsteady 3-D vortex flows. Existing methods of this type are the 
Nonlinear Discrete-Vortex methods [12-15], the Double-Panel methods [16-20], 
the Vortex-Panel method [21-22] and the Velocity-Potential -Panel method 
[23]. Flow compressibility has been accounted for by using the Prandtl- 
Glauert transformation, based on the free-stream Mach number. 

Exact integral equation formulation of the 2-0 subsonic full potential 
problem have been studied by Luu et al. [24] and Strieker [25]. However, 
these methods are not capable of dealing with flow with shocks. 

Relatively little (compared with field methods) attention has been paid, 
so far, to integral equation methods for transonic flows, certain approximate 
lEM formulations of the transonic small -perturbation problem were studied 
already in the pre-computer era, notably by Oswatitsch [26] and Spreiter 
[27], Computerized and extended versions of the approximate lEM were 
developed later by, amongst others, Norstrud [28], Crown [29] and Nixon 
[30]. It should be remarked that the approximate lEM's mentioned above, are 
all based on a special, partial-integration form of the integral equation for 
the transonic small -distrubance problem, which enables easy implementation of 
approximating assumptions on the decay of the perturbation velocity away from 
the body. The approximating assumptions on the decay of the perturbation 
velocity and the shock-fitting character of these methods are not considered 
to be competitive with FDM's. 

Later, Piers and Slooff [31] developed a lEM based on transonic small- 
distrubance theory. Their method does not contain any approximating 
assumptions and, through the introduction of artificial viscosity and 
directional bias, has a shock-capturing capability similar to that of current 
FDM's. 
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The development of the numerical methods based on the full potential 
equation has been made during the last Xmo years. Kandil and Yates [32] 
developed a lEM for steady transonic flows past delta wing, the results show 
that the method is promising and efficient. Oskam [33] developed a panel 
method for the full potential equation applied to multicomponent airfoils. 
This method is accomplished by adding a field distribution of source 
singularities to the conventional distribution of singularities over the 
boundaries of the field. At about the same time, Erickson and Strande [34] 
used Green's Third Identity as a means to extend panel method to non-linear 
potential flow, in which the concept of artificial density is employed and 
optimization technique is used to make sure the total compressibility is 
conserved. Later on, Sinclair [35] published an exact lEM for 2-D steady 
transonic flows, which is similar to Erickson's [34] method. 

The study of these integral equation methods shows a very good potential 
for lEM to challenge widely used FD methods and to replace FD methods in some 
fields of potential transonic flows, but, the development of the IE method is 
far from complete. Therefore, it is a challenge task to develop transonic 
integral equation methods. 

2. ADVANTAGES OF INTEGRAL EQUATION METHODS 

The integral equation method has several advantages over the finite 
difference method. The IE approach involves evaluation of integrals, which is 
more accurate and simpler than the FD approach in which the accuracy depends 
on the grid size. In IE approach, grid refinement and high order source- 

vorticity and compressibility modeling can be used in order to increase the 
accuracy. Moreover, the IE method automatically satisfies the far-field 
boundary conditions as 0(l/r) or 0(l/r^) and hence a small limited region 
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around the source of disturbance Is needed. In FD method, grid points are 
needed over a large region around the source of disturbance and special 
treatment is required to satisfy the far-field boundary conditions. The IE 
method is computationally inexpensive, particularly for unsteady flows, and 
does not suffer from the artificial viscosity effects as compared to FD method 
for transonic-flow-shock-capturing. 

3. FORMULATION 


The work presented here is for steady, inviscid transonic flows around 
two-dimensional configurations. 

(a) Differential Governing Equation and B.C.'s 

The governing equations of the two-dimensional, steady, inviscid 
compressible flow are given by 
Conservation Of Mass 


(p* «•**) * + 

X X 


(P *2*^ 2* ■ ® 


Conservation of Energy 


*2 ^ Y -1 *2 ^ 


*2 


Y-1 


*2 

(V 


*2 

$ *) 


U) 


( 2 ) 


Isentropic Gas equation 

p*/p: • (3) 

★ ★ ★ 
where ^ is the total velocity potential, a is the speed of sound, p is the 

density, y is the ratio of specific heats, and the subscript ® refers to the 

freestream condition. Combining equations (2) and (3) and using V^, p^ and a 

length Jl (A is the chord of airfoil) as the reference parameters, we obtain 

the dimensionless equation for the density p 
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= [1 (1 - 4 ? - 


where M is the freestream Mach number. 

CO 

The dimensionless form of equation (1) is given by 


where 


$ + $ a G 

XX 2Z 


G = -— (p $ +p_$) 

p ^X X z 


The boundary conditions are given by 


(i) Surface No-Penetration (kinematic) Condition 


V • "g = 0 


on g(r) = 0 


(ii) Kutta Condition 


A C. 


= 0 


TE 


( i i i ) I nf i ni ty Condi tion 

V $ ■»> 0 away from g 

where V is total velocity vector, g(r) = 0 is body (airfoil) 
pressure coefficient and TE refers to the trailing edge. 

(b) Integral Equation Solution 

The integral equation solution of Eq. (5) for a 
configuration is given in terms of the velocity field by 

V (x,z) = ^ * q (s) ds 

® ® u-ir + u-cr 


1 ^ .. I^\ (z-C) T - (x-c) k 

^ ^ 


ds 


ix-ir + (z-C)‘ 


(4) 

(5) 

( 6 ) 

(7) 

( 8 ) 
(9) 

surface, Cp is 
two-dimensional 


+ 


( 10 ) 
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// C) 


(x-C) i + (z-C) k 
(x-C)^ + (z-C)^ 


dCdC 


+ 



(s) 


(£2iL4— 

(x-c)^ + (2-C)^ 


ds 


where S is shock surface, q and y are surface source and vortex distribution; 
respectively. 

In Eq. (10), first integral is the contribution of the body thickness; 
second integral is the contribution of lifting or lifting and thickness or 
thickness only; third integral is the contribution due to the compressibility; 
while last integral is contribution due to shock. 

Not all of the first and second integral terms in Eq. (10) are 

necessarily included in the calculation of the velocity field. For symmetric 
flows, either first integral or second integral can be used; while for 

asymmetric flows, either second integral or both integrals should be used. 

Note that the integrand of the volume integral of Eq. (10) decreases 
rapidly with increasing distance not only because of the factor 

l/((x-C) + (z-C)*^) but also because G diminishes rapidly with increasing 

distance. Consequently, for computational purposes, the volume integral needs 
to be addressed only within the immediate vicinity of the body. 

The reader should notice that the difference between the present 
formulation and the formulations given by Sinclair [35] and by Tseng and 
Morino [36]. The present formulation is based on the velocity field in which 
the source term G contains first order derivatives of density only, and the 
normal velocity is discontinuous across the shock. The Sinclair and Tseng and 
Morino formulations are both based on the velocity potential in which the 
source term G contains first- and second-order derivatives of the velocity 
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potential and the velocity potential is continuous across the shock. The 

present formulation has two advantages over the velocity-potential 
formulation: (1) only first-order derivatives need to be calculated by finite 

differencing, and (2) one does not need to calculate derivatives of the 

velocity potential in order to detect the shock formation since the velocity 

field is calculated directly in the present formulation. 

(c) Shock-capturing and Shock-fitting Technique 

Before switching to the next section for numerical procedure it is worth 
discussing the last integral term of Eq. (10). 

For treating transonic flows, we consider two techniques - shock- 
capturing and shock-fitting techniques. Although these techniques are well 

known, we examine here their application to the integral equation formulation. 

In the shock-fitting technique, the contribution of the shock to the 
velocity field is represented by an explicit surface - integral term, the last 
term of Eq. (10), and the shock strength q^ is given by the normal velocity 

increment across it. Thus 

<5s = - <''ln - ''an) 

But in the shock-capturing technique, the last term of Eq. (10) is not 
necessarily included, since the volume integral itself implicitly includes the 
shock surface contribution [32]. 

Therefore, for shock-capturing technique Eq. (10) becomes. 


(x. 2 ) 2ui+wk = i+i-'fc q„(s) 




ds 


+ 






(s) 


ds 

ix-ir + (2-c) 


( 12 ) 
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+ -^ // G (C»C) 


(x-C) 1 + 

(x-C)^ + (2-C)^ 


d^dC 


We will use both shock-capturing and shock-fitting techniques for present 

work. 


4. SOLUTION PROCEDURE 

The whole solution procedure will be described by three parts: (1) 
discretization of equations, (ii) iterative procedure for subcritical flows, 
and (iii) treatment of critical flows. 

(i) Discretization of Equations 

Eq. (12) expresses the solution to the potential equation, Eq. (5), as the 
sum of four terms, the first three terms are the standard panel method terms 
and the last term is a volume integral in the field described above. The 

airfoil surface is divided by a number of straight panels with linear source 
and/or vortex distribution. The field is divided into a number of rectangular 
meshes except at the airfoil surface where trapezoidal elements are used. The 
value of G is taken to be constant over each mesh. Figure 1 shows the airfoil 
panels and field meshes. After discretization of equations, the integral 
equation solution, Eq. (12), becomes 


V (x,z) 





(x-C) i + 

(x-5)^ + (z-C)^ 


ds 


k 


+ 


1 

TH 


N 

z 

k=l 



(z-C) 1 - (x-C) k 

ds 

2 2 k 

(x-C) + (z-C) 


( 13 ) 
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. ^ r f G. . // ^ 

i=l k=l A.. (x-C) + (z-C) 

where N is total number of surface panels, IM x JM the total number of field 
meshes, A^-j the area of each field meshes. 

(ii) Iterative Procedure for Subcritical Flows 

The main difference between standard panel methods and integral equation 
methods (or field panel methods) is due to the volume integral term in Eq. 
(13). This term is a non-linear term and therefore unlike the standard panel 
method, the solution cannot be obtained directly and an iterative procedure is 
necessary. For simplicity, we discuss the iterative procedure for subcritical 
flows first. 


The iterative cycle is described below. 

Step 1 

A standard panel method calculation, with 6 = 0, is employed to get qg 
or Yg. 

Step 2 

Calculate 6 for each mesh by using the linear compressibility, 

2 

G = M u„, where Uy is calculated from the results of Step 1. 

Step 3 

Enforce the no-penetration condition and Kutta condition to get new qg 
or Yg* In this step, the contribution of 6 is included. 

Step 4 

Calculate 6 by using non-linear compressibility, Eq. (6). In this step, 
we first calculate $ and $ by using Eq. (13), compute density (p) by Eq. (4) 

A ^ 

and p and p are computed by central differencing for subcritical flows. 

A ^ 
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step 5 

Enforce no-penetration condition and kutta condition again to solve for 
new qg or Yg . 

Step 6 

Repeat Steps 4 and 5 until Qg or Yg and G converge. In 

each step, the pressure coefficient is computed by 

Y 

Cp = — |- {Cl + (l-V^) l'^ -1} (14) 


where V = |v| . Note that when calculating the pressure on the surface, the 
self-induced velocity must be treated carefully. 

(iii) Treatment of Critical Flows 

Several key points have been used for treating transonic flows with 
shocks: 

(1) For the critical flow with shocks, a mixed type finite difference 
must be applied to calculate and ; consistant with the mixed nuture of 
transonic flow. The local Mach number is calculated by 

M=VM^pi^ (15) 


(2) The conditions behind the shock are determined by using the Rankine- 
Hugoniot relations: 


V 


2n 


^ 

(y*1) 


(16. a) 


V 


2t “ 


V 


It 


(16. b) 
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P2 


(y+1) m 


In 


(y-1) + 2 


Pi 


(16. c) 


(3) Shock panels are introduced in order to sharpen the discontinuity in 
flow properties. The panels are added after several iterations when the 
location of the shock is fixed. The orientation of the shock panels is 
determined by Rankine-Hugoniot equation. 


p = sin 


rl.2 sin (p) sin (0) 
•- cos (p -"9) 




1/2 


(17) 


where p is the shock angle and 9 is the relative direction of the flow behind 
the shock to that ahead of the shock. The strength of shock is give by 


= - <''ln - '' 2 n> ' - II - 7->- " 

'^In 


In 


> 1.0 


(18) 


5. NUMERICAL EXAMPLES 

A scalar program has been developed to implement the solution procedure 
of the shock capturing-shock fitting technique. To prove the concept and 
verify the algorithm, the code is applied to NACA0012 airfoil at different 
Mach numbers and different angles of attack. 


CASE A Subcritical Flows 

Fig. 2 shows a comparison of the present results with experimental data 
for incompressible flow at a = 0° . A total of 140 panels is used for whole 
calculations over upper and lower airfoil surfaces. Two different modeling 
options, source panel modeling and vortex panel modeling, are applied 
respectively. The comparison shows that the vortex panel modeling is better 
than the source panel modeling. 
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Fig. 3a through Fig. 3c are for the case of a = 0, = 0.72 . Fig. 3a 

gives the comparison of the results of the two different modeling methods with 
Euler's solution, and it is again shown that the vortex panel modeling is 
better than the source panel modeling. Fig. 3b shows the effect of the size 
of computational domain, while Fig. 3c shows that the present solution 
compares well with Euler's solution. 

The next case illustrated is a lifting case with a = 2° and =0.63 . 

Fig. 4a shows the effect of the size of the computational domain and Fig. 4b 
shows the comparison of present solution with finite difference solution of 
the full potential equation. 

CASE B Critical Flows 

Fig. 5a shows the present results compared with the experimental data and 
Garabedian's finite difference solution. The computational domain used is 
2c X 1.5 c as shown in Fig. 1. The properties behind the shock are calculated 
using Rankine-Hugoniot relations as given by Eqs. (16. a) - (16. c). The 
results show that the present integral solution method is capable of capturing 
a sharp shock. The location and strength of the shock calculated here are in 
a very good agreement with the experimental data and the finite difference 
solution. Fig. 5b shows the effect of shock panels on the shock strength. 
The comparison shows that the shock panels are effective for sharpening the 
shock. Fig. 6 shows the lifting case result for a = 2°, = .75 . The 

comparison shows good agreement. 
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6. CONCLUDING REMARKS 

The standard panel method has been extended to the treatment of transonic 
flows. The results presented here show that the method is capable of 

capturing sharp shocks. The location and strength of shocks are in a very 

good agreement with the experimental data and finite difference solutions. 
Introduction of shock panels sharpens the shock with relatively coarse 
grids. Presently, the method is applied to unsteady transonic flow cases 

where the unsteady density term is treated as another volume integral term. A 
time-splitting finite difference scheme is used to step the equation in time 
while the integral solution is used to compute the spatial variation. The 
full length paper will include transonic flow cases in pitching motion. 
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cx = 0.0, with vortex “panels , 
computational domain: 2 x 1.5. 








